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SUMMARY 


This final report summarizes the research accomplished under NASA Grant NAG 2 502 
during the entire funded period which extended from February 1, 1988 to January 31, 1992. 
The total funding amounted to $173,625. The Technical Officers for this grant were Michael 
J. Green and Thomas A. Edwards of NASA Ames Research Center. 

INTRODUCTION 

The research performed during this grant represents a continuation of the work ini- 
tiated under the NASA Grants NAG 2-245 and NAS2-12861. During these grants, two 
new parabolized Navier-Stokes (PNS) codes have been developed to compute the three- 
dimensional, viscous, chemically reacting flow of air around hypersonic vehicles such as the 
National Aero-Space Plane (NASP). 

The first code (TONIC), originally developed by Prabhu et al. [1,2], solves the gas 
dynamic and species conservation equations in a fully coupled manner using an implicit, 
approximately-factored, central-difference algorithm. During the present grant, this code 
was upgraded [3] to include shock fitting and the capability of computing the flow around 
complex body shapes. The revised TONIC code was validated by computing the chemically- 
reacting (Moo=25.3) flow around a 10° half-angle cone at various angles of attack and the 
Ames All-Body model at 0° angle of attack. The results of these calculations [3] were in 
good agreement with the results from the UPS code. One of the major drawbacks of the 
TONIC code is that the central-differencing of fluxes across interior flowfield discontinuities 
tends to introduce errors into the solution in the form of local flow property oscillations. 
In order to control these oscillations, some type of artificial dissipation is required. The 
correct magnitude of this added smoothing is generally left for the user to specify through a 
trial-and-error process. 

The second code (UPS), originally developed by Lawrence et al. [4,5] for a perfect gas, 
has been extended by Tannehill et al. [6-9] to permit either perfect gas, equilibrium air, or 
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nonequilibrium air computations. The code solves the PNS equations using a finite— volume, 
upwind TVD method based on Roe’s approximate Riemann solver that has been modified to 
account for real gas effects using the approach of Grossman and Walters [10]. The dissipation 
term associated with this algorithm is sufficiently adaptive to flow conditions that, even 
when attempting to capture very strong shock waves, no additional smoothing is required. 
For nonequilibrium calculations, the code solves the fluid dynamic and species continuity 
equations in a loosely-coupled manner. The fluid medium is assumed to be a chemically 
reacting mixture of thermally perfect (but calorically imperfect) gases in thermal equilibrium. 
This code was used to calculate the hypersonic, laminar flow of chemically reacting air over 
cones at various angles of attack. In addition, the flow around the McDonnell Douglas 
generic option blended-wing-body was computed [9] and comparisons were made between 
the perfect gas, equilibrium air, and the nonequilibrium air results. 

The TONIC and UPS codes have been carefully compared during the present grant. 
During these comparison tests, [3,11] it became evident that the UPS code is much more 
robust than the TONIC code and requires less user interaction. In addition, it was noted 
that the loosely-coupled approach used in the UPS code permits different chemistry models 
to be readily inserted. As a result of the comparison tests, it was decided to discontinue the 
development of the TONIC code while continuing the development of the UPS code. 

The UPS code has since been extended in the present grant to permit internal turbulent 
flow calculations with hydrogen-air chemistry. The chemistry model contains eleven reac- 
tions and nine species and is based on the NASP model [12]. With these additions, the new 
code has the capability of computing both aerodynamic and propulsive flowfields. This ca- 
pability is required in order to successfully analyze the scramjet propulsion system employed 
on the National Aero-Space Plane. The new code has been applied to two internal flow test 
cases. The first case consists of the Burrows-Kurkov supersonic combustion experiment [13] 
in which hydrogen was injected tangentially at sonic speed through a slot in the floor of a 
test section with a Moo=2.44 vitiated airstream. In the second test case, the code was used 
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to compute a generic 3— D scramjet (M^, — 7.0) flowfield. 


GOVERNING EQUATIONS 

The PNS equations are used in the present research to model the fluid dynamics. These 
are obtained from the steady, compressible Navier-Stokes equations by neglecting stream- 
wise viscous terms and by retaining only a fraction of the streamwise pressure gradient term 
in the subsonic layer in order to eliminate ellipticity in the marching direction. The lat- 
ter is accomplished using Vigneron’s technique [14] in conjunction with the extension to 
chemically-reacting flows by Prabhu et al. [2]. The PNS equations expressed in generalized 
coordinates 3X6 gi ven by 

e* + f„ + g ( = o (1) 

where 

e - ( 7 ) e,+ ( j ) f ‘ + (^) Gi 

F = (y) (E, - e;) + (^) (Fi - f;) + ( j) (g. - g;) ( 2) 

g = (E. - e;j + (F. - f;) + (G. - g;) 

The inviscid and viscous flux vectors are given by 

T 

Ei = {pu,pu 2 +p,puv,puw,(E t + p)u} 

T 

F, = [pv,puv,pv 2 + p,pvw,(E t + p)v} 

Gi = {pw, puw, pvw, pw 2 + p,(E t + p)w} (3) 

T 

E v — {0, Tjy, T lx , UT XI 4“ "f* 

F„ = {0, Tyxt Tyy, Ty t , UTy* "f VTyy + WTy t ~ ^ 

T 

G v = {0, T zxy T Z y , T zz , tlT it + VT Z y + WT tz Q z } 

where E t = p{e -I- ^(u 2 + u 2 + to 2 )} 
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where the non-dimensional quantity #3 is 




0 * p* 
r 00 ^00 


x* Re- 


and P„„ is the multicomponent diffusion coefficient for the species a. In the present work a 
kinetic binary diffusion coefficient P is used and is assumed to be the same for all the species. 
The species continuity equation is simplified using the PNS approximation of dropping the 
unsteady term and neglecting the streamwise diffusion terms. After recasting the equation 
into generalized coordinates, the final form is obtained. 
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In addition to the above equations, the equation of state is used: 

A pT 


V = 


M 


(7) 


( 8 ) 
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where the nondimensional quantity 0 i and molecular weight of the mixture M. are given by 


0i = 


n *ik 


M 


■(Sft) 


-1 


and TZu is the unversal gas constant (8314.34 J/kmol/K). The ratio of specific heats, 7 is 
defined as 


7 = 1 + 


te) 


( 11 ) 


where e* is the sensible energy which can be expressed in terms of the species mass fractions 
(c,) and the species formation enthalpy at 0 °K (h°J by 

e* = e - J2 c ‘ h °f. 


J = 1 


In addition, the following nondimensional quantities are used 

M* 


M = 


u>: = 


M e 


c C; ' T ~ 


Ul. 


p* V 

r 00 00 


V = 


V* 

V* 


Thermodynamic and Transport Properties 


Enthalpy and specific heat 

The enthalpy and specific heat are obtained from a table lookup procedure using the 
data of Ref. [15]. Cubic spline interpolation is used to the find the property at a particular 
temperature. Since the enthalpies in Ref. [15] are referenced to 298.15 K, they are re- 
referenced to 0 °K in the following manner. For each species, the enthalpy at 0 °K is subtracted 
from the enthalpy at a particular temperature T (all referenced to 298. 15° K). This yields the 
sensible enthalpy referenced to 0 °K at the temperature T. The species formation enthalpy 
at 0°K is then added to obtain the properly referenced enthalpy. The enthalpy and frozen 
specific heat of the mixture are given by 


h * 


r , 

*r dT * 


Cl,..., Cn 


A iK 

AT- 


J=1 


( 12 ) 
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where the subscripts on the differentiation denote that the mixture composition is locally 
frozen. 

Viscosity and Thermal Conductivity 

Cubic spline interpolation is employed to obtain the species viscosity, n,, from the 
tabulated data given in Ref. [16]. The thermal conductivity of species s is computed using 
Euckeo’s semiempirical formula 


. fin- („ M\ , 5 \ 

• ~ \ n • + 4 ) 

The viscosity and thermal conductivity of the mixture are calculated using Wilke’s semiem- 
pirical mixing rule [17]. 




X./iJ 
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X, = 


c,M * 

m: 


#. = Ex. 
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Diffusion coefficient 

The binaxy Lewis number, £e, is assumed to be the same constant for all the species 
and is taken to be unity [18] for the present calculations. The kinematic diffusion coefficient 
T > * is then computed from the definition 


V* = 


K*Ce 



(14) 


Chemistry Model 

A eleven reaction/nine species hydrogen— air chemistry model is employed. The reactions 
and the corresponding forward reaction rate variables are based on the NASP model [12] 
and are given in Table 1. The forward reaction rate for the A:th reaction is expressed in the 
following expanded Arrhenius form: 

K} tk {r) = AT* n exp(-0/T*) (15) 
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Reaction 

A 

n 

0 

1 

H 

+ 

o 2 




0 

+ 

OH 

1.91E+14 

0 

8273 

2 

0 

+ 

h 2 




H 

+ 

OH 

5.06E+04 

2.67 

3166 

3 

OH 

+ 

OH 




0 

+ 

H 2 0 

1.50E+09 

1.14 

0 

4 

OH 

+ 

h 2 




H 

+ 

h 2 o 

2.16E+08 

1.51 

1726 

5 

0 

+ 

NO 




N 

+ 

o 2 

3.80E+09 

1.0 

20820 

6 

0 

+ 

n 2 




NO 

+ 

N 

1.82E+14 

0 

38370 

7 

H 

+ 

NO 




N 

+ 

OH 

1.70E+M 

0 

24560 

8 

H 

+ 

H 

+ 

M 


h 2 

+ 

M 

7.30E+17 

-1.0 

0 

9 

H 

+ 

0 

+ 

M 

- 

OH 

+ 

M 

2.60E+16 

-0.6 

0 

10 

0 

+ 

0 

4- 

M 


o 2 

+ 

M 

1.14E+17 

-1.0 

0 

11 

H 

+ 

OH 

+ 

M 

- 

H 2 0 

+ 

M 

8.62E+21 

-2.0 

0 


Table 1: Reactions and reaction rates 


In Table 1 the units for the forward reaction rates are cm 3 /mol-sec or cm 6 /mol J -sec and the 
third— body efficiencies are 2.5 for M=H 2 , 16.25 for M=H 2 0 and 1.0 for all other M. 


The above model of eleven reactions (m = 11), nine species (n = 9) and ten reactants 
(n t = 10) can be symbolically represented as: 

>1>, * = (16) 
f=l 1 

where v kl and v'h axe the stoichiometric coefficients and Ai is the chemical symbol of the 
Ith. species. Using the law of mass action, the mass production/depletion rate of the species 


r= 1 


S IS 

n n - *wn ii 

*=1 v r - 1 

The mole-mass ratios of the reactants are defined as 

Cr/M r = 1,2, . . . ,n 

r = n+l,...,n, 


7r = 


(17) 


(18) 
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where Z r> , are the third-body efficiencies for each of the species. The backward reaction rate 
required for lj* is obtained from 


K* 

is* f j i lr — 12 m 

A b t k - is* > « - 


(19) 


where K ^ t is the equilibrium constant of the fcth reaction given by 



where 7 = 82.06 x 10- 6 m 3 atm/mol-°K [19] and A n fc is the integer difference between the 
numbers of product and reactant species: 


An* = £ v' ki , - 53 v k,. 


• = 1 


J=1 


(21) 


and 


ac; = X>U - X>U 


( 22 ) 


J =1 


J=1 


The species Gibbs free energy g, are obtained from tables in Ref. [15]. In the present study 
the reactions 5,6 and 7 have been “turned off” so as to treat N 2 as an inert gas. 


Turbulence Modelling 

The algebraic turbulence model proposed by Baldwin and Lomax [20] is used in the 
present code for turbulent calculations. This model was chosen for its inherent simplicity 
and its suitability for complex flows with length scales that are not well defined. For three- 
dimensional internal corner flows, the turbulence model is modified as proposed by Hung et 
al. [21]. Using the computed eddy viscosities, the thermal conductivity and mass diffusivity 
are calculated to account for turbulent mixing. A turbulent Prandtl number of 0.9 is used 
for all the calculations in the present research. 
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NUMERICAL METHOD 


Gasdynamic Solution 

A finite-volume, upwind, TVD scheme is is used to integrate the fluid dynamic equa- 
tions. The algorithm is second-order accurate in the crossflow plane and first-order accurate 
in the streamwise marching direction. The upwind algorithm is based on Roe’s steady ap- 
proximate Riemann solver [22] which has been modified [7] for “real gas effects. The heat 
flux terms include the effects of mass diffusion. Second-order central differences are used to 
model the mass diffusion terms. Further details of the algorithm can be found in Refs. [5,9]. 

Chemistry Solution 

The species continuity equation is solved in a loosely-coupled manner using a finite- 
volume formulation. The requirement that the mass fraction of the species sum to unity 
eliminates the nth species continuity equation: 

< 23 > 

#=1 

This results in requiring only n — 1 equations to be solved. The convective terms are modeled 
using first-order upwind differences and the strong conservation— law form is retained by using 
the fluid fluxes (the coefficients of the convective terms) as known quantities from the most 
recent fluid integration step. The species production/depletion rate u> # , is treated as a source 
term and is lagged to the nth marching station for the present calculations. 

A line Gauss-Seidel procedure with successive over-relaxation (SOR) is used to solve 
each equation. A scalar tridiagonal solver is used to solve the resulting system of equations in 
an iterative manner until the residual drops below a specified tolerance level e . The residual 
is defined as 

I C*. +1 - c‘. |< £ 

where i + 1 is the current iteration level and i is the previous iteration level. 
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Fluid/Chemistry Coupling 

The coupling between the fluids and chemistry is achieved in an approximate manner. 
The fluid step is first taken with frozen chemistry to advance from the n to the n + 1 
marching station. The fluid density and velocity computed at the new station are then used 
to advance the chemistry solution to the n + 1 level. After determining the species mass 
fractions, mixture molecular weight, fluid density and internal energy at the n + 1 level, the 
new pressure, temperature, 'y, specific enthalpy and frozen specific heats are calculated. 

The temperature is obtained using the following Newton- Raphson iterative scheme 


■w*+l 



P(T* h ) - e* 

F(T* k ) 


(24) 


where 


nr k ) = i>. (wn - ^) 

Ar>) = Ec.(c;,(n-^) 

and k is the iteration level. The iterations are continued until 


(25) 

(26) 


| T k+ i _ T k |< S 

where 8 is a specified tolerance level. Once the temperature is determined, the pressure can 
be found from Eq. 10 and 7 from Eq. 11. 

The coupling between the fluids and chemistry can be enhanced through the implemen- 
tation of Newton iterations on the governing equations at each streamwise step [9]. However, 
this was found not to be necessary for the cases considered in this study. 


NUMERICAL RESULTS 

The new internal flow UPS code has been used to compute two test cases. The two test 
cases were chosen to demonstrate and validate the hydrogen— air combustion as well as the 
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three-dimensional internal flow capability of the code. The first test case is the Burrows 
Kurkov supersonic combustion experiment [13] and the second one is a three-dimensional, 
internal flow, shock induced combustion case which simulates a generic 3-D scramjet flow- 
field. 

Test Case 1: 

In the two-dimensional Burrows-Kurkov experiment [13], combustion occurs in the su- 
personic shear layer produced by the sonic injection of hydrogen into a stream of vitiated 
air. The test section consists of two nearly parallel walls with the lower wall slightly angled 
down. A schematic of the experimental setup is shown in Fig. 1. The freestream conditions 
for the hydrogen jet and the vitiated air are given in Table 2 and the wall temperature was 
held constant at 298° K. 


Freestream 

H 2 jet 

Vitiated airstream 

Mach number 

1.0 

2.44 

Temperature, °K 

254.0 

1270 

Pressure, atm 

1.0 

1.0 

H a mass fraction 

1.0 

0.0 

H 2 O mass fraction 

0.0 

0.256 

Oa mass fraction 

0.0 

0.258 

N 2 mass fraction 

0.0 

0.486 


Table 2: Freestream conditions for Burrows-Kurkov experiment 

For all the calculations, a grid consisting of 101 grid points in the normal direction was 
used. The grid was clustered near the lower wall in order to properly resolve the shear layer. 
The first point off the wall was placed at l.Ox 10~ 8 m. The Baldwin-Lomax turbulence model 

was employed to simulate turbulent mixing. 

Two computations were performed with the first one being a pure mixing case. For 
both computations, freestream startup conditions were assumed at the x=0 plane. In the 
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mixing case only, the freestream temperature was set to 1150°K to match the experiment [13] 
and all of the 0 2 in the vitiated air region was replaced by N 2 so that no combustion takes 
place. The species mole fraction profiles at the exit plane (x=35.6 cm) are compared with 
the experimental results in Fig. 2. The computed results are in excellent agreement with 
the experimental results. The second calculation used the flow conditions listed above which 
allow supersonic combustion to occur. Ignition, based on the mass fraction of OH species, 
was found to occur at about 15cm. The species mole fraction profiles at the exit plane 
are compared with the experimental results in Fig. 3. The flame strength denoted by the 
peak in the H 2 0 profile and the wall mole fraction values agree well with the experimental 
predictions but the species profiles are shifted closer to the lower wall. A series of grid 
refinement studies were performed and no appreciable change in the behavior of the profiles 
was found. This issue is currently being investigated further. The total temperature profiles 
at the exit station are compared in Fig. 4. The computed results, including the peak total 
temperature location and magnitude, compare well with the experimental data. 

Test Case 2: 

The second test case consists of a three-dimensional duct with a 15 degree compression 
ramp. Air, pre-mixed with hydrogen, enters the duct with a freestream Mach number of 
7. Combustion occurs as a result of the shock emanating from the compression ramp. The 
schematic of the 3-D duct is shown in Fig. 5. The side walls have been removed for the 
sake of clarity. The freestream flow conditions are given in Table 3. The flow is assumed 
to be turbulent and a constant wall temperature of 500°K is used. A grid consisting of 
61x61 points at each marching station was clustered at all four walls to properly resolve the 
turbulent boundary layer. Due to the presence of the strong shock, smaller streamwise step 
sizes were taken in the vicinity of the compression corner. 

The contours in the centerline, streamwise plane for Mach number and H 2 0 mole frac- 
tion are shown in Figs. 6 and 7, respectively. As expected the strong compression shock 
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Mach number 

7.0 

Reynolds number, / m 

1. 013x10® 

Temperature, °K 

1200.0 

Hj mass fraction 

0.03207 

Oa mass fraction 

0.25447 

Nj mass fraction 

0.71346 


Table 3: Freestream conditions for 3-D case 

emanating from the corner induces combustion. The profiles of pressure, temperature and 
H 2 0 mole fraction at the centerline of the exit plane (x=3cm) are shown in Figs. 8, 9 and 10, 
respectively. Once again it is clearly seen that the compression due to the shock increases 
the pressure and temperature and initiates the reactions to produce water vapor. 

CONCLUDING REMARKS 

The three-dimensional UPS code has been extended in the present grant to solve internal 
turbulent flows with hydrogen-air chemistry. The code now has the capability to compute 
external or internal flows with either perfect gas, equilibrium air, nonequilibrium air, or 
nonequilibrium hydrogen-air chemistry. As a consequence, this code can now be used to 
compute the integrated aerodynamic/propulsive flowfields of hypersonic vehicles such as the 

NASP. 
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Injection 


Figure 1: Schematic of the Burrows-Kurkov experimental setup 
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Figure 4: Total temperature profiles at the exit plane (x=35.6cm; combustion case) 
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Figure 5: Schematic of the 3-D duct 
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